# ################################################################
# ##例子5-3:控制变量方法计算欧式看涨期权价格
# T=1; r=0.04; sigma=0.3; X=5; S=8
# d1=(log(S/X)+(r+0.5*sigma^2)*T)/(sigma*sqrt(T))
# d2=d1-sigma*sqrt(T)
# V_exact=S*pnorm(d1)-exp(-r*T)*X*pnorm(d2) ##BS公式解
# 
# ##标准的蒙特卡洛方法
# N=1000; Z=rnorm(N)
# S_T=S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*Z) ##模拟到期日的股票价格
# V_M=exp(-r*T)*(S_T-X)*(S_T>X) ##计算这个股票价格对应的欧式看涨期权价格
# V_mc=mean(V_M) ##样本均值作为欧式看涨期权价格的估计值
# 
# ##控制变量蒙特卡洛方法
# rand_M=rnorm(N)
# S_T_con=S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*rand_M) ##模拟股票价格
# option_M_con=exp(-r*T)*(S_T_con-X)*(S_T_con>X) ##期权价格样本
# control_v_M=S_T_con-exp(r*T)*S ##控制变量样本
# cro_cov=cov(option_M_con,S_T_con) 
# alpha=cro_cov/var(S_T_con) ##最优系数
# V_mc_con_M=option_M_con-alpha*control_v_M ##使用控制变量方法的期权价格样本
# V_mc_con=mean(V_mc_con_M) ##均值作为估计

# ################################################################
# ##例子5-4:对偶方法计算欧式看涨期权价格
# T=2; r=0.04; sigma=0.3; X=50; S=70
# d1=(log(S/X)+(r+0.5*sigma^2)*T)/(sigma*sqrt(T))
# d2=d1-sigma*sqrt(T)
# V_exact=S*pnorm(d1)-exp(-r*T)*X*pnorm(d2) ##BS公式解
# print('BS公式计算期权价格：')
# print(V_exact)
# 
# ##标准的蒙特卡洛方法
# N=10000; Z=rnorm(N)
# S_T=S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*Z)
# V_M=exp(-r*T)*(S_T-X)*(S_T>X)
# V_mc=mean(V_M)
# print('标准蒙特卡洛方法计算期权价格：')
# print(V_mc)
# 
# ##对偶蒙特卡洛方法
# N=1000; Z1=rnorm(N)
# Z2=-Z1
# S_T_dual_1=S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*Z1)
# S_T_dual_2=S*exp((r-0.5*sigma^2)*T+sigma*sqrt(T)*Z2)
# V_M_dual_1=exp(-r*T)*(S_T_dual_1-X)*(S_T_dual_1>X)
# V_M_dual_2=exp(-r*T)*(S_T_dual_2-X)*(S_T_dual_2>X)
# V_M_dual=(V_M_dual_1+V_M_dual_2)/2
# V_mc_dual=mean(V_M_dual)
# print('对偶蒙特卡洛方法计算期权价格：')
# print(V_mc_dual)

################################################################
##习题5-2:对偶方法计算欧式看涨期权价格
##设股票价格服从几何均值回复过程 dS(t)/S(t)=a(b-S(t))dt+sdW(t)
##设 a=0.5, b=60, sigma=s=0.3
##设到期日T=2, 无风险利率r=0.04, 敲定价格X=50, 当前股价S=70

T=2; r=0.04; X=50; S=70
a=0.5; b=60; sigma=0.3
print('这个随机微分方程的公式解是什么？')

##对偶蒙特卡洛方法
N=1000 ##生成N条样本路径
M=100 ##时间区间分成100等分
dt=T/M ##每段小区间的时间长度

##产生所有路径的所有小区间上的随机增量对应的随机数
Z=matrix(rnorm(N*M),nrow=N,ncol=M,byrow=T) 

##准备存放N条样本路径的每条路径上的M+1的节点的股票价格
St=matrix(rep(0,N*(M+1)),nrow=N,ncol=M+1,byrow=T) 
##准备存放对偶路径的股票价格
Sd=matrix(rep(0,N*(M+1)),nrow=N,ncol=M+1,byrow=T) 

##每条路径的初始价格都是S=70
St[,1]=S 
Sd[,1]=S 

##准备存放N条样本路径的到期日股票价格
S_T_dual_1=rep(0,N) 
S_T_dual_2=rep(0,N) 

##使用离散格式：S(k+1)-S(k)=a(b-S(k))S(k)dt+ sigma*S(k)*sqrt(dt)*Z(k)
for (n in 1:N){
  for (k in 1:M){
    St[n,k+1]=St[n,k]+a*(b-St[n,k])*St[n,k]*dt + sigma*St[n,k]*sqrt(dt)*Z[n,k]
    Sd[n,k+1]=Sd[n,k]+a*(b-Sd[n,k])*Sd[n,k]*dt - sigma*Sd[n,k]*sqrt(dt)*Z[n,k]
  }
}

##保存N条路径的到期日股票价格
S_T_dual_1=St[,M+1]
S_T_dual_2=Sd[,M+1]

##计算N次模拟的期权价格
V_M_dual_1=exp(-r*T)*(S_T_dual_1-X)*(S_T_dual_1>X)
V_M_dual_2=exp(-r*T)*(S_T_dual_2-X)*(S_T_dual_2>X)

##标准蒙特卡洛方法估计值
V_mc=mean(V_M_dual_1)
print('标准蒙特卡洛方法计算期权价格：')
print(V_mc)

##对偶蒙特卡洛方法估计值
V_M_dual=(V_M_dual_1+V_M_dual_2)/2
V_mc_dual=mean(V_M_dual)
print('对偶蒙特卡洛方法计算期权价格：')
print(V_mc_dual)

V_max=max(max(V_M_dual_1),max(V_M_dual))
V_min=min(min(V_M_dual_1),min(V_M_dual))
par(mfrow=c(2,1))
hist(V_M_dual_1,breaks=20,xlim=c(V_min,V_max))
hist(V_M_dual,breaks=20,xlim=c(V_min,V_max))
par(mfrow=(c(1,1)))

